# coding=UTF-8

from numpy import *
def denoise(im,U_init,tolerance=0.1,tau=0.125,tv_weight=100):
    """ An implementation of the Rudin-Osher-Fatemi (ROF) denoising model
    using the numerical procedure presented in eq (11) A. Chambolle (2005).
    Input: noisy input image (grayscale), initial guess for U, weight of
    the TV-regularizing term, steplength, tolerance for stop criterion.
    Output: denoised and detextured image, texture residual. """
    """
    使用Chambolle论文中公式11的TVD去噪的数值解法
    输入：有噪声的灰度图像，初始U，TV调整项(regularizing term)的权重,步长
    输出：去噪及去纹理的图像，纹理残留
    """
    m,n = im.shape #图像的尺寸
    #初始化
    U = U_init
    Px = im #x-component to the dual field
    Py = im #y-component of the dual field
    error = 1
    while (error > tolerance):
        Uold = U
        # 计算U梯度
        GradUx = roll(U,-1,axis=1)-U # U的梯度的x分量
        GradUy = roll(U,-1,axis=0)-U # U的梯度的y分量
        #更新对偶变量（update the dual varible）
        PxNew = Px + (tau/tv_weight)*GradUx
        PyNew = Py + (tau/tv_weight)*GradUy
        NormNew = maximum(1,sqrt(PxNew**2+PyNew**2))
        Px = PxNew/NormNew # update of x-component (dual)
        Py = PyNew/NormNew # update of y-component (dual)
        #更新原始变量（update the primal variable）
        RxPx = roll(Px,1,axis=1) # right x-translation of x-component
        RyPy = roll(Py,1,axis=0) # right y-translation of y-component
        DivP = (Px-RxPx)+(Py-RyPy) # divergence of the dual field.
        U = im + tv_weight*DivP # update of the primal variable
        # update of error
        error = linalg.norm(U-Uold)/sqrt(n*m);
    return U,im-U # denoised image and texture residual

from numpy import *
from numpy import random
from scipy.ndimage import filters

#创建合成的有噪声的图像
im=zeros((500,500))
im[100:400,100:400]=128
im[200:300,200:300]=255
im=im+30*random.standard_normal((500,500))  #添加标准正态分布的噪声

U,T=denoise(im,im)
G=filters.gaussian_filter(im,10)

